1 Libraries and what not

## Base RF model
set.seed(42)

library(ggplot2)  # for autoplot() generic
library(dplyr)
library(sf)
library(stars)
library(tmap)

2 Data preprocessing

2.1 Glaciers

dat_sf <- st_read("./data/glacier_clim.shp")
## Reading layer `glacier_clim' from data source `/Users/u0784726/Dropbox/Data/devtools/glacier_mb/data/glacier_clim.shp' using driver `ESRI Shapefile'
## Simple feature collection with 95086 features and 31 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 67.47845 ymin: 27.4918 xmax: 103.8915 ymax: 45.35089
## Geodetic CRS:  GCS_unknown

2.2 Basemap

Load and crop basemap

basemap <- read_stars("~/Dropbox/Data/summer/natural_earth/HYP_50M_SR_W/HYP_50M_SR_W.tif")

basemap2 <- st_crop(basemap, st_bbox(dat_sf))
basemap2 <- st_rgb(basemap2)

3 HARR Data

projstr <- "+proj=lcc +lon_0=83.0 +lat_0=32.0 +lat_1=32.0 +lat_2=38.0 +datum=WGS84 +no_defs"

3.1 t2

t2m_2008 <- read_stars("~/Dropbox/Data/summer/harr/output/t2_2008_amjjas.nc")
st_crs(t2m_2008) <- projstr

t2m_2018 <- read_stars("~/Dropbox/Data/summer/harr/output/t2_2018_amjjas.nc")
st_crs(t2m_2018) <- projstr

Make delta

t2m_d <- t2m_2018 - t2m_2008
t2m_trend <- read_stars("~/Dropbox/Data/summer/harr/output/t2_b_ann.nc")
st_crs(t2m_trend) <- projstr

3.2 tp

tp_2008 <- read_stars("~/Dropbox/Data/summer/harr/output/prcp_2008_ann.nc")
st_crs(tp_2008) <- projstr

tp_2018 <- read_stars("~/Dropbox/Data/summer/harr/output/prcp_2018_ann.nc")
st_crs(tp_2018) <- projstr

Make delta

tp_d <- tp_2018 - tp_2008

Seasonality

tp_seas <- read_stars("~/Dropbox/Data/summer/harr/output/prcp_2018_seas.nc")
st_crs(tp_seas) <- projstr

4 Maps

4.1 t2m

x <- c(t2m_2008, t2m_2018)
names(x) <- "2m Air Temp (AMJJAS)"
m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(x) +
  tm_raster(palette = "-inferno",
            alpha = 0.75, 
            style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

Plot \(\delta\)

names(t2m_d) <- "2m Air Temp (AMJJAS) (change 2018 - 2008)"
m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(t2m_d) +
  tm_raster(alpha = 0.75, 
            style = "cont") +
  tm_layout(legend.outside = TRUE)
print(m1)

names(t2m_d) <- "2m Air Temp (AMJJAS) (change 2018 - 2008)"
m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(t2m_d) +
  tm_raster(alpha = 0.75, 
            style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

Plot \(\delta\) with mass balance (continuous scale)

m1 <- tm_shape(t2m_d) +
  tm_raster(alpha = 0.75, 
            style = "cont") +
  tm_shape(dat_sf) + 
    tm_symbols(col = "mb_mwea", 
             size = 0.25, 
             alpha = 0.75, 
             style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

4.2 t2m trend

names(t2m_trend) <- "2m Air Temp Trend (slope 1991-2020)"
m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(t2m_trend) +
  tm_raster(palette = "-PuOr",
            alpha = 0.75, 
            style = "cont", midpoint = 0) +
  tm_layout(legend.outside = TRUE)
print(m1)
## stars_proxy object shown at 1093 by 537 cells.

4.3 tp

x <- c(tp_2008, tp_2018)
names(x) <- "Total PPT (Annual)"
m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(x) +
  tm_raster(palette = "YlGnBu",
            alpha = 0.75, 
            style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

Plot \(\delta\)

names(tp_d) <- "Total PPT (Annual) (change 2018 - 2008)"
m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(tp_d) +
  tm_raster(alpha = 0.75, 
            style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

Plot \(\delta\) with mass balance

m1 <- tm_shape(tp_d) +
  tm_raster(alpha = 0.75, 
            style = "quantile") +
  tm_shape(dat_sf) + 
    tm_symbols(col = "mb_mwea", 
             size = 0.25, 
             alpha = 0.75, 
             style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

4.4 tp seasonality

names(tp_seas) <- "PPT seasonality (JJA / ANN)"
m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(tp_seas) +
  tm_raster(palette = "-PuOr",
            alpha = 0.75, 
            style = "quantile", midpoint = 0) +
  tm_layout(legend.outside = TRUE)
print(m1)
## stars_proxy object shown at 1093 by 537 cells.